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Abstract 



Operators on probability distributions can be expressed as operators on the asso- 
. dated moment sequences, and so correspond to operators on integer sequences. Thus, 

! there is an opportunity to apply each theory to the other. Moreover, probability models 

' can be sources of integer sequences, both classical and new, as we show by consider- 

■ ing the classical M/G/1 single-server queueing model. We identify moment sequences 

that are integer sequences. We establish connections between the M/M/1 busy period 
O ' distribution and the Catalan and Schroeder numbers. 

o ■ 

1 Introduction 

X. 

^ , Given the well established link between integer sequences and combinatorics, it is surprising 

that there are so few connections between integer sequences and queueing theory, because 
two prominent scholars, John Riordan (1902-1988) and Lajos Takacs (1924- ), were active 
in both combinatorics and queueing theory [22, 23, 24, 28, 29, 30]. Of course, these authors 
saw connections; e.g., the index of Riordan's queueing book [23] contains entries for Bell 
polynomials, binomial moments, Catalan numbers, cumulant generating function, Lagrange 
expansion, and rooted trees. 

In the OEIS [26], a search on "Riordan" yields 1221 entries, while a search on "queueing" 
produces 2. We think Riordan would want less disparity. One exception is the entry submit- 
ted by A. Harel (A122525) making connection to the Erlang delay formula associated with 
the classical M/M/s queue. A search on "queue" produces 15 entries, but these primarily 
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focus on combinatorial problems [11, 27] or a queue as a tool in computer science rather 
than queueing theory [12]; queueing theory concerns stochastic process models describing 
congestion, e.g., the probability distribution of customer waiting times [23, 28]. 

The purpose of this paper is to point out connections between the theories of probability 
and integer sequences, and to exhibit some integer sequences that arise in queueing theory. 
As a branch of probability theory, queueing theory has exploited the classical analytical ap- 
proach to probability theory using transforms [23, 28], whose series expansions involve the 
integer moments of the probability distributions and close relatives. The entire probability 
distribution is characterized by the sequence of moments in great generality [10, 17]. Thus 
we were motivated to introduce an operational calculus for manipulating probability distri- 
butions on the positive halfline by either manipulating the associated Laplace transforms or 
the associated moment sequences [7]; a quick overview is provided by [7, Tables 1-3]. A more 
recent paper in the same spirit is [18]. 

The operational calculus for probability distributions on the positive halfline in the frame- 
work of moment sequences is closely related to operators commonly used to analyze integer 
sequences. Since many operators on integer sequences can be applied to moment sequences 
arising in probability theory, both when they are integers and when they are not, there is 
an opportunity for experts on integer sequences to contribute to probability theory through 
moment sequences. (The connection is also discussed in [13].) The probability connection 
also provides concrete models where integer sequences arise. 

Here is how the present paper is organized. We start in §2 by reviewing moment sequences 
in probability theory and relating operators on probability distributions to operators on 
sequences. Then in §3 we review the classical M/G/1 single-server queueing model and 
identify integer sequences arising in that context. We consider the M/M/1 busy period 
distribution, its stationary excess distribution and the equilibrium time to emptiness. In that 
context we identify random quantities whose probability density functions have the Catalan 
and Large Schroeder numbers as moments. Finally, we consider the moment sequence for 
the M/G/1 steady-state waiting time. We give the proofs of Theorems 5 and 11 in §4. 

2 Moment Sequences in Probability 

Let Z be a nonnegative random variable with cumulative distribution function (cdf) F and 
probability density function (pdf) /, i.e.. 




where = means "defined as." Let f{s) be the Laplace transform (LT) of / (and thus Z) and 
let 0(x) be the associated moment generating function (mgf) of /, i.e.. 
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where s in (1) is understood to be a complex number with positive real part, while x in (1) 
is a positive real number. The LT is always well defined; we assume that there exists x* > 
such that < oo for x < x*. 

We obtain sequences by considering series expansions for the mgf 0; in particular, we can 
write 

oo „ oo 

n=0 ■ n=0 

where of course we must have /i„ = m„/n!, n > 0. From the probability perspective, 
the object of primary interest is the cdf F, but the moments m„ provide a useful partial 
characterization. The moment sequence can be calculated and employed to derive other 
quantities of interest in probability models [1, 7, 16]. There is a developing operational 
calculus for manipulating probability distributions via their moment sequences [7]. 

From the sequence perspective, we may instead regard the sequences {m„ : n > 1} and 
: n > 1} as the objects of primary interest. The two perspectives meet with the mgf 0. 
From the sequence perspective, (f){x) arises as the generating function (gf) of the sequence 
{fin '■ n > 1}, and we speak of {/in} as being its coefficients, while 0(x) arises as the 
exponential generating function (egf) of the sequence {m^ : n >1}. Of course, {m„ : n > 1} 
is always an integer sequence whenever {/i^ : n > 1} is, but not necessarily conversely. We 
will be giving examples where both are integer sequences. 

A simple example from probability theory is the exponential distribution with mean m, 
specified by 

F{t) = 1 - e"*/'", f{t) = (l/m)e-*/'", t > 0, and 0(x) = (1 - mx)-\ 

which has associated sequences 

m„ = n!m" and fin = f^^i n>Q. 

It is immediate that {/i„} and are both elementary integer sequences whenever m is an 

integer. Since the mean (first moment) is probabilistically only a scale parameter, depending 
on how the measurement units are defined, it is natural to follow the convention that the 
first moment is 1; here that gives mi = 1. Then we obtain the fundamental integer sequences 
m„ = n\ and /i„ = 1, ?7, > 1. 

As a consequence of the relations outlined above, we see that results about probabil- 
ity distributions can be translated into results about integer sequences, provided that the 
sequence {fin} or {m„} is indeed an integer sequence. Similarly, results about integer se- 
quences can be translated into results about probability distributions, provided that the egf 
or gf of the integer sequence is indeed the mgf of a bonafide pdf. 

We first observe that there is a natural probabilistic setting in which, not only is {m„} a 
moment sequence, but so is the associated sequence {fin}- That occurs in spectral represen- 
tations, where one pdf serves as a mixing pdf for another. In particular, suppose that the 
pdf / can be represented as a continuous mixture of exponential pdf's via 

f{t)= r y^'e~'lyg{y)dy, t > 0, (3) 

J ri 
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in which case we call g the mixing pdfioi f [4]. Append a superscript / to {m„} and {/in} 
to denote dependence upon /. We observe that the sequence {/i„} is itself the moment 
sequence of the associated mixing pdf g. Hence we call {/in} the mixing moments of /. (We 
omit the elementary proofs in this section.) 

Proposition 1. {mixing moments) For pdf's f and g related via (3), = /i^, n > 1. 

A canonical operation in probability theory is convolution. If two independent nonneg- 
ative random variables Zi and Z2 with pdf's fz^ and are added, then the sum Zi + Z2 
has a pdf that is the convolution of the two component pdf's, i.e.. 



Jo 

One reason transforms are so frequently used in probability is that they convert convolution 
into a simple product; i.e., the associated mgf's are related by 



Because of the assumed stochastic independence, the moments are related by the binomial 



Thus, if {m^^} and {m^^} are integer sequences, then so is {m^^~^^^}. For example, by 
above, that occurs when Zi and Z2 have exponential or deterministic distributions with 
integer means. 

With this background, we can interpret [7, Tables 1-3], which show various operators 
mapping one probability distribution into another. There are four columns: The first column 
contains the name and notation for the operator; the second column shows how the operator 
acts on LT's; the third column shows how it acts on pdf's; and the fourth column shows how 
it acts on moment sequences. From the perspective of integer sequences, the fourth column 
shows how it acts on the coefficients of the egf; we could then add a fifth column showing 
how it acts on the corresponding coefficients of the gf. Those familiar with integer sequences 
might want to translate the LT into the mgf by replacing s with —x, and then interpret that 
mgf as the egf of the given k^^ moment. 

In the rest of this section we highlight a few striking connections between operators 
on probability distributions, as in [7], and operators on integer sequences. First, a standard 
operator on integer sequences is the simple shift to the left, e.g., converting 1, 1, 2, 6, 22, 90, . . . 
to 1, 2, 6, 22, 90, . . .. We now show that, probabilistically, the simple shift applied to the 
coefficients /in corresponds to constructing the stationary-excess cdf of a cdf on the positive 
halfline having mean 1. 

Given a nonnegative real-valued random variable Z with cdf F having finite moments 
nik, k > 1, let Ze be a random variable with the associated stationary- excess cdf Fg (a.k.a. 
the equilibrium excess or stationary residual-life cdf), defined by 




theorem. 




F,{t) = P{Z, < t) 



mi Jo 



[ {1-F{u))du, t > 



(4) 
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The stationary-excess cdf frequently arises in renewal theory; see of [25, Examples 7.16,7.17, 
7.23, 7.24] and [31]; it appears in [7, Table 2]. (A search on "renewal theory" in the OEIS 
gives two unrelated entries.) For us, the important fact is that the random variable has 
moments 

Hence, the transformation from a cdf of a nonnegative random variable to its associated 
stationary-excess cdf produces a simple shift on the gf coefficients. 

Proposition 2. [the stationary- excess operator) Let F be the cdf of a nonnegative random 
variable with mean 1, mgf (p in (1) and associated sequence of mixing moments {^n} in (2) 
[coefficients of (f){x) when it is regarded as a gf). The associated stationary- excess cdf in 
(4) has mgf (f)f.{x) = {(f>{x) — l)/x. The mixing moments of 4>e{x) {coefficients of (j)e when it 
is regarded as a gf) are fie,k = /^fc+i > k > 1. 

A similar relationship holds for the stationary-lifetime operator, mapping a pdf / into 
the associated pdf 

= t>0. (6) 

mi 

The stationary-lifetime pdf also frequently arises in renewal theory, e.g., [25, §7.7], and also 
appears in [7, Table 2]. For us, the important fact is that the moments are related by 

ms,k = k>l, 
nil 

just like (5) without the + 1 in the denominator. Hence, the transformation from a pdf 
of a nonnegative random variable to its associated stationary-lifetime pdf produces a simple 
shift on the moments m„ (the egf coefficients). 

Proposition 3. {the stationary- lifetime operator) Let f be the pdf of a nonnegative random 
variable with mean 1, mgf (p in (1) and associated sequence of moments {m„} in (2). The 
associated stationary-lifetime pdf fs in (6) has mgf (j)s{x) = (j)'{x). The coefficients of (ps 
when it is regarded as an egf are rris.k = fnk+i, k > 1. 

We now turn to another basic probability operator, which is conveniently related to a 
continued fraction representation of the mgf [8, 9, 13, 15]. If / is the LT of a pdf /, then the 
associated exponential mixture pdf has LT 

f£M{s)^{l + sf{s))-'; (7) 

it appears in [7, Table 3]. The special case of exponential mixtures of inverse Gaussian 
(EMIG) distributions is discussed in [7, §8] and in [9]. 

The probabilistic exponential mixing operator has a simple manifestation in the continued 
fraction representation of the mgf, regarding that mgf as a gf. Starting with the LT f{s) of 
a pdf /, if we represent the associated mgf as a formal power series by 

(j){x) = /(-x) = 1 + fiiX + /i2X^ + /isX^ + fi^X^ + /isX^ + . . . , (8) 
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then the corresponding continued fraction (CF) is 



1 hix h-^x 

<P{x) = - — — — (9) 

where hi = fii, h2 = (yU2 — fiD/fii, etc. When > for all n we have an S-fraction and 
the underlying pdf / is completely monotone (CM) [8]. 

Proposition 4. {the exponential-mixture operator) Let f be the pdf of a nonnegative random 
variable with mean 1, mgf (p in (1) with associated sequence of CF coefficients {hn} in (9). 
The associated exponential-mixture pdf f^M ^i^h LT in (7) has CF coefficients 



h£M,n — hn-i, n >2; 



i.e., it produces a shift to the right. The inverse exponential mixture operator [7, (7.3), p. 
94] gives the corresponding shift to the left. 



Proof. From the transform expression in (7), the conclusion is immediate: Given the CF 
representation for f{s), it is immediate that the corresponding CF for (1 + sf{s))~^ shifts 
the coefficients one to the right; i.e., we write the CF for (1 + sf{s))~^ as 1/(1 + sf{s)), 
inserting the CF for f{s). m 



3 Queueing Examples 

3.1 The M/G/1 Model 

The M/G/1 queue is a basic model in queueing theory, usually discussed in queueing text- 
books; e.g., [25, §8.5]. There is a single server with unlimited waiting room. Customers 
arrive according to a Poisson process (the first M, for Markov) with rate A, < A < oo. 
If the system is empty, then the customer goes immediately into service; otherwise the cus- 
tomer waits in queue. The successive service times come from a sequence of independent and 
identically distributed (i.i.d.) random variables with cdf G having mean 1///, < < oo. 
We will mostly consider the easiest case, in which the service-time cdf G is exponential; then 
the model is denoted by M/M/1. 

A waiting customer enters service immediately upon service completion. Let Q{t) be 
the number of customers in the system at time t for t > 0. In the M/M/1 model, the 
stochastic process Q = {Q(t) : t > 0} is a birth-and-death stochastic process, with constant 
birth rate A and constant death rate fi. Let p = X/fi he the traffic intensity. If p < 1, 
then P{Q{t) = j\Q{0) = z) — j- (1 — p)p^ as t — )■ oo for each i and j; i.e., Q{t) converges in 
distribution to a geometric distribution on the nonnegative integers, having mean p/{l — p). 
We assume that p < 1, under which the system is said to be stable. (If p > 1, then 
P{Q{t) < j|(5(0) = z) as t oo for each i and j.) For the M/G/1 model, the steady- 
state distribution of Q{t) is characterized by the PoUaczek-Khintchine transform [23, 28]. 
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3.2 The M/M/l Busy Period 



A busy period is the time from the arrival of a customer finding an empty system until the 
system is empty again [23, §4.8]. The first passage time of Q from any state j > to j — 1 
is distributed as a busy period. Without loss of generality, we can measure time in units of 
mean service times, so that we let /i = 1. Then the model has the single parameter p (= A). 
Let Xp denote a random variable with the busy period distribution, as a function of the 
traffic intensity p. Let Bp be the cdf of Xp, i.e., Bp{t) = P{Xp < t), t > 0, and let bp be the 
associated pdf. It turns out that 

bp{t)^-^e~^'+^^'h{2t^), t>0, 
where hit) is the Bessel function of the first kind [23, (39), p. 63]. The LT of bp (and thus 

1 + P + s - v/(l + P + s)2 -4p 



Xp) is 



bp{s) = / e~%{t)dt = E[e 
Jo 



2p 

[23, (38), p. 63]. As usual, the moments can be obtained by differentiating the Laplace 
transform. The first two moments are -E[Xp] = (1 — p)~^ and -E[^p] = 2(1 — p)~^. 



3.3 The Busy-Period Moment Sequence After Scaling 

Our goal here is to obtain interesting integer sequences from the sequence of successive 
integer moments of a busy period Xp. To obtain integer sequences, we first perform a 
change of variables, introducing a = p/(l — p) (the mean steady state number in system) or, 
equivalently, p = cr/(l + a). Notice that the first two moments become 1 + a and 2(1 + cr)^. 
It turns out that the entire moment sequence becomes an integer sequence whenever a is an 
integer. Hence, from this first step, we obtain an entire sequence of integer sequences. 

To seek simple integer sequences, we further scale these moment sequences all to have 
mean (first moment) 1. We remark that this final spatial scaling also plays a role in under- 
standing the performance of the M/M/l queue as the traffic intensity p increases toward its 
critical value 1. With appropriate scaling of both time and space, the stochastic process Q 
approaches refiected Brownian motion with negative drift, while the busy period distribution 
has interesting behavior, in which both small values and large values play a role [3, 6, 32]. 

In terms of random variables, let 

l^ = ^^[i±^ = (l-p)X„ where p = and a = (10) 
l + fj 1 + a 1 — P 

Let b{x] cr) be the mgf of y^-. From above, 

b{x;(r) = E[e''^^] = l + ^^~^~^(^) ^ ^^lere ^(x) = v^l - 2{l + 2a)x + x^. (11) 

2a 

Let the associated moments of Y„ be mn{a) = E[Y^], n > 1, where mi (a) = 1 for all a > 0. 
These moments, divided by n\, are the coefficients of the series expansion of b{x; a) 

b{x;a) = } '-f — => bn{a)x'', where bn{a) = (12) 

n=0 n=0 
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The coefficients 6ri(o") (and thus also the moments m„(o")) are polynomials in a (and thus 
integers when a is an integer); the ffist few are: 60 ('^) = 1) ^i(^) = 1? ^2(0") = 1 + 0" and 
h^[a) = 1 + 3(T + 2(t2. 

More specifically, we now relate the coefficients 6n(o") in (12) to the Catalan numbers, 
denoted by C^, (A000108) starting with 1, 1, 2, 5, 14, 42, 132, 429, 1430; in particular. 



Cn = —^\] and Cn+i = y^CiCn-h n>l. 
n + l\ n I ^-^ 

The Catalan numbers can be characterized in terms of their generating function 



(13) 



cfe).5:c./.i^4i5. (14) 

fc=o 



Theorem 5. {mixing moments of the busy period pdf) For n > 1, 

k=0 ^ ^ 

where bnicr) is defined in (11) and (12), and Ck is the k^^ Catalan number in (13). In 
addition, there is a convenient recurrence relation, with bo{a) = fci(cr) = 1 and 

(n + l)bn+i{a) = {2n - 1)(1 + a)6„(a) - (n - 2)6„_i(a). (16) 

We provide a proof in §4. For cr = 1, 2, 3, the coefficient sequences are: 

{bnil)} = 1, 1, 2, 6, 22, 90, 394, (A155069) 
{bn{2)} = 1, 1, 3, 15, 93, 645, 4791, (A103210) 
{6„(3)} = 1, 1, 4, 28, 244, 2380, 24868, (y4103211) 

Sequence {&„(1)} (A155069) is a relatively recent addition to the OEIS [26], having the title 
"Expansion of (3 — x — Vl — 6x + x^)/2 in powers of x." We provide a model context. 
However, sequence {6„(1)} shifted one to the left, i.e., 1,2,6,22,90,... is (A006318), cor- 
responding to the famous "Large Schroeder numbers." The "little Schroeder numbers" in 
(A001003) are obtained by dividing A006318 by 2, i.e., the sequence 1,1,3, 11, 45, . . .. 

The moment sequences themselves, m„ (cr) = n!6„((j) are of course also integer sequences, 
but they are evidently not in the OEIS [26]. For example, the sequences with terms n!5„(l) 
and n!6„(l/2), n>0, yielding 1, 1,4,36,528, 10800 and 1, 1,3, 18, 171,2250, respectively, are 
not found. 

The busy period is the first passage time of Q from 1 to 0. Since the first passage time 
pdf from each state /c to is the fc-fold convolution of the busy period pdf [3, Theorem 3.1]), 
the moments and mixing moments of these pdf's also generate integer sequences for integer 
cr. 

Since we find regularity by multiplying {6„(cr)} by 1/2 = o-/{l + a) when a = 1, we 
are motivated to consider the associated sequences {abn+iic) / {1 + a) : n > 1} for integers 
a > 1. 
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Corollary 6. For n > 1, 

k=0 



1 + a ^ \n — k 



We have already observed that Corollary 6 gives the little Schroeder numbers for o" = 1; 
it is also discussed in [24, Problem 15, p. 168]. 

3.4 The Busy-Period Stationary-Excess Distribution 

Let Bf, denote the busy period stationary-excess cdf associated with the busy-period cdf Bp, 
defined as in (4) after applying the scaling in (10). Let be be the associated stationary-excess 
pdf. We can apply Proposition 2 to characterize the mixing moments 6e,n(o") of fee- 
Corollary 7. {mixing moments of the busy-period stationary- excess pdf) The M/M/1 busy- 
period stationary excess pdf fee has mgf 

be{x, a) = r e-\{t) dt = Y, &e,n(a)a;" = ^^""^ ~ ^ 

•^0 n=0 ^ 

2 



1-X+ v/l-2(l + 2a)x + x2 
c{ax/{l-xf) 2 



X 1 — X + '^{x) 



(17) 



where c{y) is the generating function of the Catalan numbers in (14) and \1/ is defined in 
(11). Forn > 1, 

fee,n(cr) = bn+l{a) , U > 1, (18) 

for which an expression is given in Theorem 5. The mean of be is a. 

3.5 Large Schroeder and Catalan Numbers as Moments of PDF's 

In this section we identify the pdf's whose moments are (i) the Large Schroeder numbers and 
(ii) the Catalan numbers. By Corollary 7 and our observation after Theorem 5, the sequence 
{fee,n(l)} coincides with the Large Schroeder numbers (A006318). From [3, Theorem 5.1 and 
Corollary 5.2.1] and [4, Theorem 4.1], we can obtain the spectral representation for the pdf 
fee. (We need to account for the different scaling of time used there.) That identifies the 
desired pdf and, by Proposition 1, the mixing moments. 

Corollary 8. {Large Schroeder numbers as moments) The large Schroeder numbers arise as 
the moments of the mixing pdf associated with the stationary- excess busy-period pdf be when 
a = 1. The mixing pdf for be as a function of a is 



r, ^ v/(r-y)(y-r-i) 1 

f{y\<y) = - 7, ' -<y<r, 

zany r 



where r = 1 + 2(T + 2^a{l + a) 
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We introduced the Catalan numbers before Theorem 5. The Catalan numbers are related 
to reflected Brownian motion with drift —1 and diffusion coefficient 1, denoted by {R(t) : 
t > 0}. It is the limit of the M/M/1 queue length process as p t 1 with appropriate scaling 
of time and space [2, 32]. Since E[R{t)\R{0) = 0] is nondecreasing, Hi{t) = E[R{t)\R{0) = 
0]/E[R{oo)], t > 0, is a cdf. By [2, Corollary 1.5.2], [7, §7] and [8, (8.7)], the pdf hi oi Hi has 
LT hi{s), which can be characterized as the unique fixed point of the exponential mixture 
operator, i.e., 

hi{s) = 1— . (19) 

1 + shi{s) 

Together with the spectral representation given in [4, Theorem 4.1], that implies the following 
result. Without the probability model context, the result already appears in the commentary 
on (A000108) [21]. 

Theorem 9. {Catalan numbers as moments) The generating function c{x) of the Catalan 
numbers in (13) and (14) coincides with the mgf of hi, hi{—x), associated with the LT in 
(19). As a consequence, the Catalan numbers arise as the moments of the mixing density of 
hi, 

f(y) = ^^-y < y < 4. (20) 

If we take the Laplace transform of / in (20) and replace s by —x, then we obtain the 
mgf, which is the egf c{x) of the Catalan numbers. 

Corollary 10. The egf of the Catalan numbers is 

c{x) ^ V ^ = e^^{h{2x) - Ii{2x)), 

n=0 

where Iq and Ii are Bessel functions. 



3.6 The M/M/1 Equilibrium Time to Emptiness 

From a sequence perspective, the shifting by 1 and multiplying by cr/(l + a) following 
Theorem 5 and Corollary 6 lead us to consider the mgf 

1 (y ( b(x: a) — l\ , , 

p{x; a) ^ —— + —— -^-^ . 21 

1 + cr l + cr\ X J 

The function p{x; a) in (21) turns out to be the moment generating function of the equi- 
librium time to emptiness in the M/M/1 queue, i.e., the first passage time to state by the 
stochastic process Q, assuming that Q starts at time according to its (geometric) steady- 
state distribution. Hence, with probability (1 — p) = 1/(1 + cr) the process starts at 0, so 
the first passage time is 0. As a consequence, there is an atom at with mass 1/(1 + cr). 
Conditional on starting at a positive value, the equilibrium time to emptiness coincides with 
the stationary-excess pdf be of the busy-period pdf b [5, Theorem 3]. Additional charac- 
terizations appear in [3, Theorem 3.3]. We make a connection to the exponential mixture 
operator in (7). Paralleling (12), we let Pn(c") be the coefficient of p{x] a) as a gf, i.e., writing 
p(x; cr) = J2'^=QPn{o-)x^. We provide a proof in §4. 
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Theorem 11. {equilibrium time to emptiness) The mgfp{x;a) in (21) can be expressed as 

= c((l + a)x/(l+xn ^ 2 

^^"^^ ' 1 + x l + x + vl>(x)' 

where c{y) is the generating function of the Catalan numbers in (14) and \1/ is defined in 
(11). Hence, 

Pn{cr) = ^ = — — > J \Ck(T , n > 0, 23 

1 + cr 1 + a ^-^ \n — k J 

k=o ^ ^ 

and the following recursion can be applied with po{a) = 1 and pii^cr) = cr 

(n + 2)p„+i(or) = (2n + l)(l + 2a)p„,(a)-(n-l)p„_i(a), n > 1. (24) 
In addition, b{x;a) can be expressed as an exponential mixture ofp{x;a), i.e., 

h{x-a) = - ^. (25) 

1 — xp[x] a) 

From (23) and (24), the sequence {pn{o') : n > 1} is an integer sequence for each positive 
integer a. For a = 1,2,3, the coefficient sequences are: 

{Pnil)} = 1, 1, 3, 11, 45, 197, (A001003) 
{p„(2)} = 1, 2, 10, 62, 430, 3194, (A107841) 
{p„(3)} = 1, 3, 21, 183, 1785, 18651, (A131763). 

We conclude this section by remarking that in the theory of integer sequences there also 
is a convenient invert operator, which can be expressed for LT's via 

g{s) = X{f{s)) = Excess(exp mixture(/(s)) = fsM,e{s) 

Excess I — I = - I 1 



,^ + sf{s)) sy l + sf{s)^ 
^^'^ (26) 



l + sf{s 



The following result links the three mgf's h{x] cr), he{x\ cr) and (p(x; cr); see Corollary 5.1 
of [3]. We also consider the excess of the excess, denoted by he.e [31]. 

Corollary 12. The mgf's be{x;a), p{x;a) and b{x;a) are related by 

be{x; cr) = X{p{x; cr)) = p{x; a)b{x; a) 
for the invert operator in (26) and 

be,e{x;a) =Peix; a). 
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Directly from the final expression in (17), we get a heavy-traffic limit for the scaled mgf 
with our scaling in (10), with the limit being the gf of the Catalan numbers. Let Xf,{a) be a 
random variable with mgf be{x;a); i.e., be{x;a) = Ele^-^^^"'^]. Let =^ denote convergence in 
distribution [32]. 

Corollary 13. As a cx). 



bJx/a; a) , = = cix) 

for c in (14 ), so that 



where c{x) = E[e^^]. 

A related heavy-traffic limit for the time-scaled busy-period stationary-excess pdf was 
obtained in [6, Theorem 1]. (They are consistent.) 

3.7 Continued Fractions and Hankel Transforms 

We now look at the three probability mgf's b{x; a), be{x, ; a) = {b{x; a) — l)/x and p{x; a) in 
(11), (17) and (21) from the perspective of continued fractions and Hankel transforms. From 
the previous subsections, we know that the associated three random variables are closely 
related. That is shown again through this new perspective. First, the three mgf's can be 
represented as S (Stieltjes) fractions. (In [8] we found that S fractions frequently arise in 
CF's associated with birth-and-death stochastic processes.) The CF coefficients are given in 
Table 1. 



mgf 


hi 


h2 


h3 


hi 


/l5 


he 


b{x; a) 


1 


a 


1 + a 


a 


1 + a 


a 


be{x] a) 


1 + ^7 


a 


1 + a 


a 


1 + a 


a 


p{x; a) 


a 


1 + a 


a 


1 + a 


a 


1 + a 



Table 1: The coefficients of the continued fractions representing the three M/M/1 mgf's. 
The pattern repeats in each row. 



XJa) 



X as (J — )■ 00, 



a 



For the generating functions of our M/M/1 queueing examples, the determination of the 
CF coefficients is simple because we can apply the algebraic identity 



a 7 a 7 

T+T+T+T+'" 2 



1 (v/l + 2(7 + a) + (7-a)' - 1 - (7 - «)) 



for constants a and 7; i.e., we can solve the equation 

a 7 



c 



1+1 + c 

12 



for C- 



The Hankel transform of an integer sequence provides a useful partial characterization; it 
is a many-to-one function mapping an integer sequence into another integer sequence. (For 
an example of non-uniqueness, see Corollary 14 below.) For example, the Hankel transform 
of the Catalan numbers is the sequence 1, 1, 1, . . . [13, 20]. Following [15, 17], starting from 
a sequence {cj„ : n > 0} = Wq, Wi, cj2; i^2; • • • with = 1, let the Hankel matrix M^"^ be the 
(n + 1) X (n + 1) symmetric matrix with elements M}"^ = Ci;j+j_2; < i < n, < j < n. 
(The first row contains the first n -\- 1 elements and Mn+i^n+i = ^2n- Let H2n = rfet(M(")), 
the even Hankel determinant. Let the Hankel transform of the sequence {a;„ : n > 0} above 
be the sequence {H2n '■ n > 0}; it starts with Hq = 1. 

One way to compute the Hankel transform of a sequence {un : n > 0} is to first determine 
the corresponding continued fraction from the formal power series; i.e., starting from (8), we 
determine (9) above. The coefficients hn appearing in (9) can be determined by applying 
the normalized Viskovatov algorithm, as given on [15, p. 112]. Then we apply the iteration 



see Theorem 1.4.10 on of [17, p. 23]. Another way to arrive at this result is via the even 
contraction of the CF in (9), as given in [13, (12.3), p. 270]; note that f3n in [13] is equal 
to h2n-ih2n, > 1, in our notation. Then (12.2) of [13, (12.2)] and our iteration yield the 
same result. 

Corollary 14. The Hankel transforms of the sequences {&n(o")} in Theorem 5 , {&e,n(c")} in 
Corollary 7 and {Pn{(^} in Theorem 11 are 



where z/(n) = n{n + l)/2. 

3.8 The Stationary Waiting Time in the M/G/l Queue 

The mgf w{x\ a) = E[e^'^] of the steady-state waiting time W in the M/G/l queue can be 
expressed in terms of the mgf ge{x]a) of the stationary-excess of the general service-time 
distribution using the PoUazcek-Khintchine transform (again assuming that the mean service 
time is 1) by applying the random sum operator from [7, Table 1], yielding 

w{x; a) = ^ , = — —, r (27) 

[23, §4.3] and [5]. Riordan [23, (18a), p. 49] and Takacs [28] develop a nice recursion for the 
moments, which ClS db function of a becomes 




H2n{h) 
H2n{p) 



H2n{h) = {a + aY^^\ 




(28) 
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where is the k^^ moment of the service-time cdf. Clearly, (28) can be the source of many 
integer sequences. 

To illustrate, we now consider the Catalan numbers as service-time moments, which 
is legitimate by Theorem 9. We exploit the following result, obtained by combining [2, 
Corollaries 1.3.2 and 1.5.1] with Theorem 9. 

Lemma 15. Let hi be the density of the RBM first moment function, which has the Catalan 
numbers as its mixing moments. Then 

hi^eS = hi{sY, so that Ce{x) = c{xY, (29) 

where c{x) is the generating function of the Catalan numbers. 

Theorem 16. // the service time pdf is the RBM first moment pdf hi characterized by its 
LT in (19), which has the Catalan numbers as its moments, then the waiting time mgf in 
(27) becomes 

w{x;a) = - --TT (30) 

and the moments of the stationary waiting time pdf are 



E[W-~\a)] = ^E7!r:^C',E[^"^'(^)]- (31) 

k=2 



{n-k)\ 

Hence, if a is an integer, then {E[W^{a)]} is an integer sequence. 



For the case a = 1, we find {E[iy"(l)]} = 1, 2, 18, 252, 4776, . . ., which is not in the OEIS. 
Riordan [23, (19), p. 50] shows how the moments -E[W^"(o")] can be expressed in terms of 
the multivariate Bell polynomials. 

Let Wnic) = E[W"'{a)]/n\ be the associated mixing moments. From (31), we obtain a 
recursion for Wn(o"). 

Corollary 17. A recursion for the mixing moments w„((j) defined above is 

n 

w„_i((t) = (y^CkWn-k{(y), n>l. (32) 

k=2 

For cr = 1, we get {tf„(l)} = 1, 2, 9, 42, 199, . . ., which also is not yet in OEIS. 



3.9 General Birth-and-Death Processes 

The most relevant previous work connecting queueing theory to integer sequences and pro- 
viding a suitable framework for generalization evidently is the previous work connecting 
birth-and-death processes to continued fractions; see [8, 18] and references therein; contin- 
ued fractions are known to be intimately connected to integer sequences. We have already 
mentioned that the stochastic process Q = {Q(t) : t > 0}, representing the number of 
customers in an M/M/1 queueing model at time t, is a birth-and-death stochastic process. 
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Indeed, it is a special birth-and-death process with constant birth rate A and constant death 
rate /i. More generally, these rates are functions of the state [25, §6.3]; is the rate up and 
Hk is the rate down when Q{t) = k. The general birth-and-death process represents a more 
general queueing model. 

For a general birth-and-death process X = {X[t) : t > 0}, an important quantity is the 
probability 

Po,o(t) = P{X{t) = 0\X{0) = 0), t> 0; (33) 

see [3, §11], [18, §3.1] and [19]. From [18, (4.1), p. 771], the LT of Po,o{t), Po,q{s), has a 
representation as an S fraction 

^o,o[s) — —7- — -7- — 

S+1+S+1+S + 
_ Z XqZ (liZ XiZ fX2Z ^^^-^ 



1+ 1+ 1+ 1+ 1- 



where z = 1/s. (The last line follows from [8, (1.5)] using the sequence {c„ : n > 0} with 
C2n = 1 and C2n+i = z, n > 1.) We can find the corresponding power series 

-Po,o(s) = z{l -piz+p2Z^ -Psz^ + ■■■); 

see [8, (3.9), (3.10), p. 397]. The sequence {1,^1,^25^37 • • •} Hiay be an integer sequence. 

For the special case of the M/M/1 queue, with our scaling we have Xi = a and /Xj+i = 1+a 
for all i, i > 0, so that 

- / , f I az (1 + a)z az (1 + a)z 

Then from Table 1 we have the relation 

Po,o{s) = zp{z), (35) 
where p{s) is the Laplace transform of the equilibrium time to emptiness. 



4 Proofs 

In this concluding section we provide proofs for Theorems 5 and 11. 



4.1 Proof of Theorem 5. 

We give two proofs. The first is direct; the second applies [14]. 

Direct proof of (15). From [29, pp. 232-233], after a change of scale and notation, we 
have 
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(This equation is also given by [24, (21), p. 151] after we make the identification that 
Hn{z) = Pnio') = cr6„4.i(cr)/(l + 0").) On the right side, move (1 + cr) inside the sum and 
identify the coefficients of a'', obtaining 

/ 1 \ / /'n + k + l\/'n-l\ /'n + k\/n-l\\_/n + k\/2k\/' 1 \ 
\n+l) VV ^ + 1 )\ k k )\k-l) ) ~ {n-k){k ) [ITT) ' ' 

Direct proof of (16). We follow the technique of [23, p. 107] to establish a three-term 
recursion. By the successive differentiation with respect to x in (11), we find that 

^{xfb"{x; (t) = (1 + (t)(1 + 2a- x)b'{x; a) + (1 + a)%{x; a). 

In this equation, make the following substitutions: 

oo 

b{x; a) = ^6„(a)a;", 

n=0 

oo 

b'{x; a) = ^(n + l)bn+ii(j)x'^, 

n=0 

oo 

b"{x; a) = 5^(n + l){n + 2)6„+2(o-)x". 

n=0 

Then collect and equate the coefficients of x^~^ and the result follows. ■ 



Application of [14]. The idea in this second proof is to directly make connection to the 
Catalan numbers via their generating function c in (14). Starting from (11), we apply (17), 
(21) and the proof of (22) in Theorem 11 below to conclude that 

We then apply [14, Proposition 15, p. 12] (with change of variables y = ax and y = {1 + a)x, 
respectively) to immediately deduce the conclusions 

&e,n(a) = X^(^;^J;)c,a^ and P.(a) = ^ ^ ^) (-1)"^'=C,(1 + a)^ (37) 

A;=0 ^ ^ fc=0 ^ 

given in Corollary 7 and Theorem 11. We then can apply (18) and (17) to get the conclusion 
in Theorem 5 for b from b^- By this line of reasoning, we establish Theorems 5 and 11 and 
Corollaries 6 and 7 all at once, and have a new perspective on their relationship. 
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4.2 Proof of Theorem 11. 

Proof of (22). From the representation of p{x; a) in (21), we have 

p{x]a) 



1 + 0- 1 + a \1- x + "^{x) 
1 a /l-x-'^(x) 



1 + 0- 1 + o- V 2ax 
l + x-^(x) 2 



2(l + a)a; l + x + '^{x)' 



On the other hand, by (14), 

c((l + a)a;/(l + x)2) 



l + x 1 + a; + (1 + x) v^l - 4(1 + a)x/{l + x^ 

2 

~ l + x + -^{x)' 

Proof of (23). Combine (15) and (21). ■ 

Proof of (24). Combine (16) and (21). ■ 

Proof of (25). From (11), 

1 + 2(T - X - ^Cx) 2(1 + a) 



b{x; a) 



2a 1 + 2o- - X + \[^(x) 

1 1 



1 - X ([1 + X - ^(x)]/[2(l + o-)x]) 1 - xp(x; a) ' 
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